# Author: Marco Portmann (marco.portmann@unifr.ch)
# Please ask permission before using this skript in your projects
# Last change: 12.11.2015
#
#



mmp <- function(regressors=NA, pretty.regressors=NA,
                allmodels,
                EmptyCols=0,
                decimals=4,
                filename="clipboard",
                # out of order
                LaTeX=F,
                CoefOnly = F,

                n.rgroup=0, rgroup=NA,
                bottom.matter=NA, models.names=NA,
                landscape=F,  caption, caption.location="bottom",
                label="tab:", append=F, factors = NA, coefaspercent = F) {

# to do: rgroups=NA
#        models von summary(obj) auf obj ändern
#        ... Kompatibilität zu glm, survreg, robcov etc.
#        Intercept (Intercept) ?

require(Hmisc)
require(lmtest)


# summary statistics: variable names
# first col.: use this names in the function argument list
# second col.: text which appears in the table
# third col.:  text which appears in a LaTeX table
# ... could add columns with long/short names

mod.info <-
matrix(
c(
"sigma"                     ,  "Sigma"             ,  "$\\sigma$"         ,
"DF"                        ,  "DF"                ,  "DF"                ,
"adj.r.squared"             ,  "Adj. R2$"          ,  "Adj. $R^2$"        ,
"f-statistics"              ,  "F-stat."           ,  "F-stat."           ,
#"DF"                       ,  "DF"                ,  "DF"                ,
"Breusch-P."                ,  "Breusch-P."        ,  "Breusch-P."        ,
"N"                         ,  "n. Obs."           ,  "n. Obs."           ,
"r.squared"                 ,  "R2"                ,  "$R^2$"             ,
#"PseudoR2"                  ,  "PseudoR2"          ,  "Pseudo $R^2$"     ,
"LogLik"                    ,  "Log-Likelihood"    ,  "Log-Likelihood"    ,
"Brier"                     ,  "Brier"             ,  "Brier"
  ), byrow=T, ncol=3)

if (LaTeX==T) info.col <- 3 else info.col <- 2



# use fixed number of decimals (given by var "decimals"), -->> ###use -e0x for small numbers
  printnumber <- function(num, ispercent = F){
      if (ispercent==T) num <- num * 100
      spstring <- paste("%.",decimals,"f", sep="")
      res <- sprintf(spstring, num)
      if (is.na(num)==F){
        if (abs(num)<=10^(-(decimals-1))) {
           spstring <- paste("%.",1,"e", sep="")
           res <- sprintf(spstring, num)
        }
      }else{res <- "NA"}

      if (ispercent==T) res  <- paste(res, "%", sep="")

      return(res)
   }



# significance stars
  GetStars <- function(t) {
    t <- abs(t)
    n <- -1 + as.numeric(cut(t,breaks=c(-0.01,-qnorm(c(0.05, 0.025, 0.005)),Inf)))
    if (n == 0) {
      return("")
    } else {
      return(paste("", paste(rep("*", n), sep="", collapse=""), "", sep=""))
    }
  }

# significance stars
  stars.tstat <- function(t, df) {
    t <- abs(t)
    n <- -1 + as.numeric(cut(t, breaks=c(-0.01,-qt(c(0.05, 0.025, 0.005),df),Inf)))
    #print(n)
    if (n == 0) {
      return("")
    } else {
      return(paste("", paste(rep("*", n), sep="", collapse=""), "", sep=""))
    }
  }




#.........................................................
#
# Extract model information from different classes
#

 models <- list(rep(NA, length(allmodels)))
for (j in 1:length(allmodels)) {
      m <- allmodels[[j]]

      if (is.null(m)) models[[j]] <- list(rep(NULL, 10))
       
      if (class(m)[1] == "lm"){

        mm <- as.matrix(
                            cbind(
                                  summary(m)$coefficients[,1],   # Coef
                                  summary(m)$coefficients[,2],   # Std. Error
                                  summary(m)$coefficients[,3],   # t value
                                  summary(m)$df[2]               # Stars
                                 ))

        loc <- grep("(Intercept)",names(m$coefficients))
        if (length(loc)>0) names(m$coefficients)[loc] <- "Intercept"

        rownames(mm) <- names(m$coefficients)


        models[[j]] <-
                      list(
                            mm,
                            c(
                                  summary(m)$sigma   ,             #"Sigma"
                                  summary(m)$df[2]   ,             #"DF"
                                  summary(m)$adj.r.squared,        #"Adj. R2"
                                  summary(m)$fstatistic[1],        #"f-statistic"
                                  bptest(summary(m))[4],           #"Breusch-P."
                                  length(residuals(summary(m))),   #"n. Obs."
                                  summary(m)$r.squared,            #r.squared"
                                  "",                              #"Log-Likelihood"
                                  ""                               #"Brier"
                             ))

       }

      if (class(m)[1] == "boot"){

        mm <- as.matrix(
                            cbind(
                                  m$t0,   # Coef
                                  apply(m$t, 2, sd),   # Std. Error
                                  m$t0/apply(m$t, 2, sd),  # Normal
                                  -1                                 # Stars
                                 ))

        loc <- grep("(Intercept)",names(m$t0))
        if (length(loc)>0) names(m$t0)[loc] <- "Intercept"

        rownames(mm) <- names(m$t0)


        models[[j]] <-
                      list(
                            mm,
                            c(
                                  NA   ,             #"Sigma"
                                  NA   ,             #"DF"
                                  NA,        #"Adj. R2"
                                  NA,        #"f-statistic"
                                  NA,           #"Breusch-P."
                                  nrow(m$data)                      , # "n. Obs."
                                  NA,            #r.squared"
                                  "",                              #"Log-Likelihood"
                                  ""                               #"Brier"
                             ))

       }

      if (class(m)[1] == "ols"){

        mm <- as.matrix(
                            cbind(
                                  m$coefficients,                    # Coef
                                  sqrt(diag(m$var)),                 # Std. Error
                                  m$coefficients/sqrt(diag(m$var)),  # Normal
                                  -1                                 # Stars
                                 ))
        rownames(mm) <- names(m$coefficients)


        models[[j]] <-
                      list(
                            mm,
                            c(
                                  ""                                , # "Sigma"
                                  ""                                , # "DF"
                                  ""                                , # "Adj. R2"
                                  ""                                , # "f-statistic"
                                  ""                                , # "Breusch-P."
                                  m$stats[[1]]                      , # "n. Obs."
                                  m$stats[[4]]                      , # "r.squared"
                                  ""                                , # "Log-Likelihood"
                                  ""                                  # Brier

                             ))
       }


      if (class(m)[1] == "tobit"){

        mm <- as.matrix(
                            cbind(
                                  m$coefficients,                    # Coef
                                  sqrt(diag(m$var))[1:length(m$coefficients)],                 # Std. Error
                                  m$coefficients/sqrt(diag(m$var))[1:length(m$coefficients)],  # Normal
                                  -1                                 # Stars
                                 ))
        rownames(mm) <- names(m$coefficients)


        models[[j]] <-
                      list(
                            mm,
                            c(
                                  ""                                , # "Sigma"
                                  ""                                , # "DF"
                                  ""                                , # "Adj. R2"
                                  ""                                , # "f-statistic"
                                  ""                                , # "Breusch-P."
                                  ""                                , # "n. Obs."
                                  ""                                , # "r.squared"
                                  m$loglik                          , # "Log-Likelihood"
                                  ""                                  # Brier

                             ))
       }



      if (class(m)[1] == "lrm" & (class(m)[2] == "rms" | class(m)[2] == "Design") & class(m)[3] == "glm"){

        mm <- as.matrix(
                            cbind(
                                  m$coefficients,                    # Coef
                                  sqrt(diag(m$var)),                 # Std. Error
                                  m$coefficients/sqrt(diag(m$var)),  # Normal
                                  -1                                 # Stars
                                 ))
        rownames(mm) <- names(m$coefficients)


        models[[j]] <-
                      list(
                            mm,
                            c(
                                  ""                                , # "Sigma"
                                  m$stats[[4]]                      , # "DF"
                                  ""                                , # "Adj. R2"
                                  ""                                , # "f-statistic"
                                  ""                                , # "Breusch-P."
                                  m$stats[[1]]                      , # "n. Obs."
                                  m$stats[[10]]                     , # "r.squared"
                                  m$stats[[3]]                      , # "Log-Likelihood"
                                  m$stats[[11]]                       # Brier

                             ))
       }

#############

      if (class(m)[1] == "mer" ){

        mm <- as.matrix(
                            cbind(
                                  fixef(m),                          # Coef
                                  sqrt(diag(vcov(m))),               # Std. Error
                                  fixef(m)/sqrt(diag(vcov(m))),      # Normal
                                  -1                                 # Stars
                                 ))
        rownames(mm) <- names(fixef(m))


        models[[j]] <-
                      list(
                            mm,
                            c(
                                  ""                                , # "Sigma"
                                  ""                                , # "DF"
                                  ""                                , # "Adj. R2"
                                  ""                                , # "f-statistic"
                                  ""                                , # "Breusch-P."
                                  "xxx???"                          , # "n. Obs."
                                  ""                                , # "r.squared"
                                  logLik(m)[1]                      , # "Log-Likelihood"
                                  ""                                  # Brier

                             ))
       }






###########


      if (length(class(m))>1){
      if (class(m)[1]=="data.frame" & class(m)[2] == "DiscreteEffects" ){

        mm <- as.matrix(
                            cbind(
                                  m[,1],                  # Coef
                                  m[,2],                  # Std. Error
                                  m[,1]/m[,2],            # Normal
                                  -1                      # Stars
                                 ))
        rownames(mm) <- sub('^DE', "", rownames(m) )


        models[[j]] <-
                      list(
                            mm,
                            c(
                                  ""                              , # "Sigma"
                                  ""                              , # "DF"
                                  ""                              , # "Adj. R2"
                                  ""                              , # "f-statistic"
                                  ""                              , # "Breusch-P."
                                  ""                              , # "n. Obs."
                                  ""                              , # "r.squared"
                                  ""                              , # "Log-Likelihood"
                                  ""                                # Brier

                             ))
       }}






      if (length(class(m))>1){
      if (class(m)[1]=="data.frame" & class(m)[2] == "CoefTest" ){

        mm <- as.matrix(
                            cbind(
                                  m$Difference,           # Coef
                                  m$SE,                   # Std. Error
                                  m$Normal,               # Normal
                                  -1                      # Stars
                                 ))
        rownames(mm) <- sub('^DE', "", rownames(m) )


        models[[j]] <-
                      list(
                            mm,
                            c(
                                  ""                              , # "Sigma"
                                  ""                              , # "DF"
                                  ""                              , # "Adj. R2"
                                  ""                              , # "f-statistic"
                                  ""                              , # "Breusch-P."
                                  ""                              , # "n. Obs."
                                  ""                              , # "r.squared"
                                  ""                              , # "Log-Likelihood"
                                  ""                                # Brier

                             ))
       }}




      if (length(class(m))>1){
      if ((class(m)[1]=="numeric" | class(m)[1]=="matrix") & class(m)[2] == "SummaryStat" ){

        mm <- as.matrix(
                            cbind(
                                  m,           # Coef
                                  -1,                   # Std. Error
                                  -1,               # Normal
                                  -2                      # Stars
                                 ))
        rownames(mm) <- names(m)
        if (is.null(names(m))) rownames(mm) <- rownames(m)

        models[[j]] <-
                      list(
                            mm,
                            c(
                                  ""                              , # "Sigma"
                                  ""                              , # "DF"
                                  ""                              , # "Adj. R2"
                                  ""                              , # "f-statistic"
                                  ""                              , # "Breusch-P."
                                  ""                              , # "n. Obs."
                                  ""                              , # "r.squared"
                                  ""                              , # "Log-Likelihood"
                                  ""                                # Brier

                             ))
       }}




      if (length(class(m))==1){
      if (class(m)[1]=="data.frame"){

        mm <- as.matrix(
                            cbind(
                                  m[,1],                  # Coef
                                  NA,                     # Std. Error
                                  NA,                     # Normal
                                  -2                      # Stars
                                 ))
        rownames(mm) <- rownames(m)


        models[[j]] <-
                      list(
                            mm,
                            c(
                                  ""                              , # "Sigma"
                                  ""                              , # "DF"
                                  ""                              , # "Adj. R2"
                                  ""                              , # "f-statistic"
                                  ""                              , # "Breusch-P."
                                  ""                              , # "n. Obs."
                                  ""                              , # "r.squared"
                                  ""                              , # "Log-Likelihood"
                                  ""                                # Brier

                             ))
       }}



      if (class(m)[1] == "negbin"){
      
        # Adaption for covariance matrix provided by the user, e.g. for 
        # cluster robust standard erros.
        # Package checks whether the $var exists otherwise it darws on
        # vcov().
        
        if(is.null(m$var)){usevar <- vcov(m)}else{usevar <- m$var}
      
      
        mm <- as.matrix(
                            cbind(
                                  m$coefficients,                    # Coef
                                  sqrt(diag(usevar)),                # Std. Error
                                  m$coefficients/sqrt(diag(usevar)), # Normal
                                  -1                                 # Stars
                                 ))

        loc <- grep("(Intercept)",names(m$coefficients))
        if (length(loc)>0) names(m$coefficients)[loc] <- "Intercept"

       rownames(mm) <- names(m$coefficients)


        models[[j]] <-
                      list(
                            mm,
                            c(
                                  ""                                , # "Sigma"
                                  ""                                , # "DF"
                                  ""                                , # "Adj. R2"
                                  ""                                , # "f-statistic"
                                  ""                                , # "Breusch-P."
                                  length(m$y)                       , # "n. Obs."
                                  m$stats[[4]]                      , # "r.squared"
                                  ""                                , # "Log-Likelihood"
                                  ""                                  # Brier

                             ))
       }



  }

#
#.........................................................

#xxx <<- models

# create list with all variables if regressors is empty

# begin with the most complete model
for(i in  order(rapply(models, length, classes="matrix"), decreasing = T)){
  if (exists("all.regressors")==F){
    all.regressors <- rownames(models[[i]][[1]])
  }else{
        if (length(rownames(models[[i]][[1]])[(rownames(models[[i]][[1]]) %in% all.regressors)==F]))
        all.regressors <-c(all.regressors, rownames(models[[i]][[1]])[(rownames(models[[i]][[1]]) %in% all.regressors)==F])
       }
}

if (length(regressors) != length(pretty.regressors) & (length(pretty.regressors)==1 & mode(pretty.regressors)=="logical") ==F ){
    warning("Length of pretty.regressors does not match length of regressors")
    pretty.regressors <- NA
}


if (length(regressors)==1)
  if (is.na(regressors[1])==T){
    regressors <- all.regressors
    pretty.regressors <- NA
    }

# use pretty.regressors or regressors if pretty.regresors is not provided
if (length(pretty.regressors)==1 & mode(pretty.regressors)=="logical") {
  pretty.regressors <- regressors
}


factors.all <- data.frame(FactorName = as.character(), FactorFound = as.character(), stringsAsFactors = F)
# Delete factors and add fixed effects
if (length(factors)>0 & is.na(factors[1]) == F){
 
 for (i in 1:length(factors)){
   factorlist <- unique(c(
                          grep(paste("^", factors[i], "=.", sep = "" ), all.regressors),
                          grep(paste("^factor\\(", factors[i], "\\).", sep = "" ), all.regressors)
                         ) 
                       )
   if (length(factorlist) > 0)   
   {
    factors.all <- rbind(factors.all, data.frame(FactorName = factors[i],  FactorFound = all.regressors[factorlist]))  
    dellist <- regressors %in% all.regressors[factorlist] == F 
    regressors <- regressors[dellist]
    pretty.regressors <- pretty.regressors[dellist]
    
   }
   
   #if (length(factorlist)>0 & sum(regressors == factors[i]) == 0)
   if (sum(regressors == factors[i]) == 0)   
   {regressors <- c(regressors,  factors[i])
      pretty.regressors <- c( pretty.regressors ,  factors[i])
     }  
 }
}
# print(factors.all)

# Warn if some variables are not printed
unused <- all.regressors[all.regressors %in% c(regressors, as.character(factors.all$FactorFound))==F]
if (length(unused[is.na(unused)==F])>0)
  warning("Unused variables: ", list(unused[is.na(unused)==F]))


 mdls <<- models



################################################################################
#
#
# no need to write interaction terms as \\* - this is done here:
regressorsX <<- regressors
if (length(grep("\\\\\\*", regressors))>0){
 stop("''\\*'' not allowed in regressors")
}

regressors <- gsub("\\*",  "\\\\*", regressors)
regressorsX2 <<- regressors
#
################################################################################



# build matrix

length.bm <- length(bottom.matter)

if (length(bottom.matter)==1) if (is.na(bottom.matter)) length.bm <- 0

modelcols <- 1+ EmptyCols
modelrows <- 1
numbers <- matrix("", nrow = length(regressors)
                           + length.bm,

                      ncol = modelcols * length(allmodels))
if (CoefOnly==T){
  numbers <-  as.data.frame(
                            matrix(as.numeric(0), nrow = length(regressors)
                                   + length.bm,

                            ncol = modelcols * length(allmodels))
                           )


}



##  colnames(numbers) <- models.names
if (is.data.frame(numbers)==F){
  rownames(numbers) <- rep("", nrow(numbers))
  colnames(numbers) <- rep("", ncol(numbers))
}


  baserow <- 1
  for (i in 1:length(regressors)) {
    basecol <- 1

    # lm: Intercept is written in brackets!
    #if (regressors[i] == "(Intercept)") regressors[i] <- "Intercept"

    if (regressors[i] == "Intercept")
     regex <- paste(regressors[i], sep="") else{
       regex <- paste("^",regressors[i],"$", sep="")
#       regex <- sub("\\*", "\\\\*", regex)
       }
    rownames(numbers)[baserow] <- pretty.regressors[i]


  
  for (j in 1:length(allmodels)) {

      m <- models[[j]][[1]]

      
      #print(regex)
      # print model name (if supplied by user)
      if (length(models.names)>=j) colnames(numbers)[basecol] <- models.names[j]
      #
      #print(regex)
      #print( rownames(m))
      #print(grep(regex, rownames(m)))
      #print(j)
      # Check whether coefficient exits and add coefficient to the table
      if (any(locations <- grep(regex, rownames(m)))){
        if (length(locations) > 1) {
          cat("Regex ", regex, " has multiple matches at model ", j, "\n")
          return(NULL)
	      }


        if (m[locations, 4]!= -2){
          if (m[locations, 4]== -1)
          pstars <- GetStars(m[locations, 3])
            else pstars <- stars.tstat(m[locations, 3], m[locations, 4])

          if (CoefOnly==F){

          thiscoefaspercent <- coefaspercent
          if(length(coefaspercent) == 1 | length(coefaspercent)<j){
           thiscoefaspercent <- coefaspercent
          }else{
            thiscoefaspercent <- coefaspercent[j]
          }

          numbers[baserow, basecol] <- paste(
                                          printnumber(m[locations, 1], thiscoefaspercent),
                                          pstars,
                                          "\n(" ,
                                          printnumber(m[locations, 2]),
                                          ")"
                                          , sep="")
          }else{
          numbers[baserow, basecol] <- (m[locations, 1])
          }
        }else{

          if (CoefOnly == F) numbers[baserow, basecol] <-   printnumber(m[locations, 1])
          if (CoefOnly == T) numbers[baserow, basecol] <-   as.numeric(m[locations, 1])


        }
      }

      if ((length(locations) == 1) &  length(grep(regex, factors))>0){
          stop( paste("Factor name ", regex, " coincides with variable name.", sep = ""))}

      if ((length(locations) == 0) &  length(grep(regex, factors))>1){
          stop( paste("Found multiple instances of factor ", regex, ".", sep = ""))}


      if ((length(locations) == 0) &  length(grep(regex, factors))==1){
       
       #if (length( grep(paste("^", regressors[i], sep = ""), rownames(m))   ) > 0)
       # numbers[baserow, basecol] <- "yes" else numbers[baserow, basecol] <- "no"
       if (nrow(subset(factors.all, FactorName == regressors[i] & FactorFound %in%  rownames(m)))>0)    
         numbers[baserow, basecol] <- "yes" else numbers[baserow, basecol] <- "no"
              
      }

    basecol <- basecol + modelcols
    }
    baserow <- baserow + modelrows

  }

# bottom.matter

  if (length.bm>0)
  if (is.na(bottom.matter)[1]==F)
  for (i in 1:length(bottom.matter)) {
    basecol <- 1


    regex3 <- paste("^",bottom.matter[i],"$", sep="")



    locations <- grep(regex3, mod.info[,1])
    if (length(locations) == 0){
          cat(bottom.matter[i], "is undefined!", "\n")
          return(NULL)

	     }
    #print(i)
    rownames(numbers)[baserow] <- mod.info[locations, info.col]
    #print('-')
    for (j in 1:length(allmodels)) {
      if (length(models[[j]])>1)
      {
        m <- models[[j]][[2]]
        
        tmp <- as.numeric(m[locations])
        if (is.na(tmp)==F){
         numbers[baserow, basecol] <- printnumber(tmp)
        }
      }
    basecol <- basecol + modelcols
    }
    baserow <- baserow + 1

  }
    
  print(numbers)

  if (filename=="clipboard")
    write.table(numbers,"clipboard", sep="\t", col.names = NA , row.names = TRUE)
#  else{
#
#  }


  return(numbers)
}




mkSummaryStat <- function(mat, fnctn){
  res <- apply(mat, 2, fnctn)
  class(res) <- c(class(res), "SummaryStat")
  return(res)
}

mkPseudoStat <- function(mat){
  mat <- t(mat)
  class(mat) <- c(class(mat), "SummaryStat")
  return(mat)
}




formatnumber <- function(numb, tstat, decimals = 4){

  printnumber <- function(num){
      spstring <- paste("%.",decimals,"f", sep="")
      res <- sprintf(spstring, num)
      if (is.na(num)==F){
        if (abs(num)<=10^(-(decimals-1))) {
           spstring <- paste("%.",1,"e", sep="")
           res <- sprintf(spstring, num)
        }
      }else{res <- "NA"}
      return(res)
  }



# significance stars
  GetStars <- function(t) {
    t <- abs(t)
    n <- -1 + as.numeric(cut(t,breaks=c(-0.01,-qnorm(c(0.05, 0.025, 0.005)),Inf)))
    if (n == 0) {
      return("")
    } else {
      return(paste("", paste(rep("*", n), sep="", collapse=""), "", sep=""))
    }
  }



return(paste(
             printnumber(numb),
             GetStars(tstat),
             "\n(" ,
             printnumber(tstat),
              ")"
              , sep="")
      )
}